rm(list = ls())
load("dataBJPOLS.RData")

inst.1 <- "judgeiv_hd * as.factor(race4) * as.factor(median_income_group)"
endo.1 <- "pti * as.factor(race4)  * as.factor(median_income_group)"
outc.1 <- "vote_post"

time.controls <- "as.factor(court_time1) + as.factor(court_time2) + as.factor(court_dow) + as.factor(severity)"
demo.controls <- "age + I(age^2) +  as.factor(race4) + female + vote_pre + as.factor(noteli) + regis_before"
case.controls <- "as.factor(any_drug) +  as.factor(any_weapon) +  as.factor(any_prop) + as.factor(any_prior_case)"

form.1 <- formula(paste(outc.1, "~", endo.1, "+" , time.controls, "|", inst.1, "+", time.controls))
form.2 <- formula(paste(outc.1, "~", endo.1, "+" , time.controls, "+", demo.controls, "|", inst.1, "+", time.controls, "+", demo.controls))
form.3 <- formula(paste(outc.1, "~", endo.1, "+" , time.controls, "+", demo.controls, "+", case.controls, "|", inst.1, "+", time.controls, "+", demo.controls, "+", case.controls))

resultsB <- resultsH <- resultsW <- resultsIB <- resultsIU <- resultsHIB <- resultsHIU <- resultsWIB <- resultsWIU <- list()
boot.run <- T

if(boot.run) {
    set.seed(1231)
    for(i in 1:500) {
        out <- last.cases[, this_id_fa[sample.int(.N, .N, TRUE)], by = c("judge_cat", "court_time1")]
        last.cases.boot <- last.cases[last.cases$this_id_fa %in% out$V1, ]
        
        m1a1 <- ivreg(form.1, data = last.cases.boot)
        m1a2 <- ivreg(form.2, data = last.cases.boot)
        m1a3 <- ivreg(form.3, data = last.cases.boot)
        
        resultsB[[i]] <- c(m1a1$coefficients["pti"], 
                           m1a2$coefficients["pti"], 
                           m1a3$coefficients["pti"])
        
        resultsW[[i]] <- c(m1a1$coefficients["pti:as.factor(race4)W"], 
                           m1a2$coefficients["pti:as.factor(race4)W"], 
                           m1a3$coefficients["pti:as.factor(race4)W"])
        
        resultsH[[i]] <- c(m1a1$coefficients["pti:as.factor(race4)H"], 
                           m1a2$coefficients["pti:as.factor(race4)H"], 
                           m1a3$coefficients["pti:as.factor(race4)H"])
        
        resultsIB[[i]] <- c(m1a1$coefficients["pti:as.factor(median_income_group)Below"], 
                           m1a2$coefficients["pti:as.factor(median_income_group)Below"], 
                           m1a3$coefficients["pti:as.factor(median_income_group)Below"])

        resultsIU[[i]] <- c(m1a1$coefficients["pti:as.factor(median_income_group)Unavailable"], 
                            m1a2$coefficients["pti:as.factor(median_income_group)Unavailable"], 
                            m1a3$coefficients["pti:as.factor(median_income_group)Unavailable"])
        
        resultsHIB[[i]] <- c(m1a1$coefficients["pti:as.factor(race4)H:as.factor(median_income_group)Below"], 
                            m1a2$coefficients["pti:as.factor(race4)H:as.factor(median_income_group)Below"], 
                            m1a3$coefficients["pti:as.factor(race4)H:as.factor(median_income_group)Below"])

        resultsHIU[[i]] <- c(m1a1$coefficients["pti:as.factor(race4)H:as.factor(median_income_group)Unavailable"], 
                            m1a2$coefficients["pti:as.factor(race4)H:as.factor(median_income_group)Unavailable"], 
                            m1a3$coefficients["pti:as.factor(race4)H:as.factor(median_income_group)Unavailable"])
        
        resultsWIB[[i]] <- c(m1a1$coefficients["pti:as.factor(race4)W:as.factor(median_income_group)Below"], 
                            m1a2$coefficients["pti:as.factor(race4)W:as.factor(median_income_group)Below"], 
                            m1a3$coefficients["pti:as.factor(race4)W:as.factor(median_income_group)Below"])
        
        resultsWIU[[i]] <- c(m1a1$coefficients["pti:as.factor(race4)W:as.factor(median_income_group)Unavailable"], 
                            m1a2$coefficients["pti:as.factor(race4)W:as.factor(median_income_group)Unavailable"], 
                            m1a3$coefficients["pti:as.factor(race4)W:as.factor(median_income_group)Unavailable"])
    }
    
    SEsB <- round(apply(do.call('rbind', resultsB), 2, sd), 3)
    SEsW <- round(apply(do.call('rbind', resultsW), 2, sd), 3)
    SEsH <- round(apply(do.call('rbind', resultsH), 2, sd), 3)
    SEsIB <- round(apply(do.call('rbind', resultsIB), 2, sd), 3)
    SEsIU <- round(apply(do.call('rbind', resultsIU), 2, sd), 3)
    SEsHIB <- round(apply(do.call('rbind', resultsHIB), 2, sd), 3)
    SEsHIU <- round(apply(do.call('rbind', resultsHIU), 2, sd), 3)
    SEsWIB <- round(apply(do.call('rbind', resultsWIB), 2, sd), 3)
    SEsWIU <- round(apply(do.call('rbind', resultsWIU), 2, sd), 3)

    ## Black x Income
    SEsBIA1 <- round(apply(do.call('rbind', resultsB), 2, sd), 3)
    SEsBIB1 <- round(apply(do.call('rbind', resultsB) + do.call('rbind', resultsIB), 2, sd), 3)
    SEsBIU1 <- round(apply(do.call('rbind', resultsB) + do.call('rbind', resultsIU), 2, sd), 3)
    
    ## Hispa x Income
    SEsHIA1 <- round(apply(do.call('rbind', resultsB) + do.call('rbind', resultsH), 2, sd), 3)
    SEsHIB1 <- round(apply(do.call('rbind', resultsB) + do.call('rbind', resultsIB) + do.call('rbind', resultsH) + do.call('rbind', resultsHIB), 2, sd), 3)
    SEsHIU1 <- round(apply(do.call('rbind', resultsB) + do.call('rbind', resultsIU) + do.call('rbind', resultsH) + do.call('rbind', resultsHIU), 2, sd), 3)
    
    ## White x Income
    SEsWIA1 <- round(apply(do.call('rbind', resultsB) + do.call('rbind', resultsW), 2, sd), 3)
    SEsWIB1 <- round(apply(do.call('rbind', resultsB) + do.call('rbind', resultsIB) + do.call('rbind', resultsW) + do.call('rbind', resultsWIB), 2, sd), 3)
    SEsWIU1 <- round(apply(do.call('rbind', resultsB) + do.call('rbind', resultsIU) + do.call('rbind', resultsW) + do.call('rbind', resultsWIU), 2, sd), 3)
}

m1a1 <- ivreg(form.1, data = last.cases)
m1a2 <- ivreg(form.2, data = last.cases)
m1a3 <- ivreg(form.3, data = last.cases)

m1a1_d <- summary(m1a1, vcov = sandwich, diagnostic = T)
m1a2_d <- summary(m1a2, vcov = sandwich, diagnostic = T)
m1a3_d <- summary(m1a3, vcov = sandwich, diagnostic = T)

res1_pti <- coefficients(m1a1)[grep("pti", names(coefficients(m1a1)))]
res2_pti <- coefficients(m1a2)[grep("pti", names(coefficients(m1a2)))]
res3_pti <- coefficients(m1a3)[grep("pti", names(coefficients(m1a3)))]

## Black Defendant: Below
c1_bb <- c(res3_pti["pti"] + res3_pti["pti:as.factor(median_income_group)Below"], SEsBIB1[3])

## Hispa Defendant: Below
c1_hb <- c(res3_pti["pti"] + res3_pti["pti:as.factor(median_income_group)Below"] + 
               res3_pti["pti:as.factor(race4)H"] + res3_pti["pti:as.factor(race4)H:as.factor(median_income_group)Below"], SEsHIB1[3])

## White Defendant: Below
c1_wb <- c(res3_pti["pti"] + res3_pti["pti:as.factor(median_income_group)Below"] + 
               res3_pti["pti:as.factor(race4)W"] + res3_pti["pti:as.factor(race4)W:as.factor(median_income_group)Below"], SEsWIB1[3])


## Black Defendant: Above
c1_ba <- c(res3_pti["pti"], SEsBIA1[3])

## Hispa Defendant: Above
c1_ha <- c(res3_pti["pti"] + res3_pti["pti:as.factor(race4)H"], SEsHIA1[3])

## White Defendant: Above
c1_wa <- c(res3_pti["pti"] + res3_pti["pti:as.factor(race4)W"], SEsWIA1[3])


m1 <- data.frame(rbind(c1_bb,
                       c1_hb,
                       c1_wb))

colnames(m1) <- c("estimate", "std.error")
m1$term <- c("t1", "t2", "t3")
m1_df <- data.frame(as.matrix(m1)) 
m1_df$estimate <- as.numeric(as.character(m1_df$estimate))
m1_df$std.error <- as.numeric(as.character(m1_df$std.error))

g2 <- {dwplot(m1_df, 
              dot_args = list(size = 2.5), whisker_args = list(size = 0.75),  
              vline = geom_vline(xintercept = 0, colour = "grey60", linetype = 2)) %>% 
    relabel_predictors(c(t1 = "Black", t2 = "Hispanic", t3 = "White")) +
    theme_classic() + xlab("The Effect of Pretrial Incarceration on Turnout") + ylab("Defendant's race") + xlim(-.25, .25)+
    ggtitle("Below Median Income") +
    theme(plot.title = element_text(face="bold"), legend.position = "none",
          axis.text.x  = element_text(size = 16),
          axis.text.y  = element_text(size = 16),
          text = element_text(size = 16))} 

m1 <- data.frame(rbind(c1_ba,
                       c1_ha,
                       c1_wa))

colnames(m1) <- c("estimate", "std.error")
m1$term <- c("t1", "t2", "t3")
m1_df <- data.frame(as.matrix(m1)) 
m1_df$estimate <- as.numeric(as.character(m1_df$estimate))
m1_df$std.error <- as.numeric(as.character(m1_df$std.error))

g3 <- {dwplot(m1_df, 
              dot_args = list(size = 2.5), whisker_args = list(size = 0.75),  
              vline = geom_vline(xintercept = 0, colour = "grey60", linetype = 2)) %>% 
        relabel_predictors(c(t1 = "Black", t2 = "Hispanic", t3 = "White")) +
        theme_classic() + xlab("The Effect of Pretrial Incarceration on Turnout") + ylab("Defendant's race") + xlim(-.25, .25)+
        ggtitle("Above Median Income") +
        theme(plot.title = element_text(face="bold"), legend.position = "none",
              axis.text.x  = element_text(size = 16),
              axis.text.y  = element_text(size = 16),
              text = element_text(size = 16))} 

pdf(file = "./Figures/Figure_03_BJPOLS.pdf", width = 12, height = 6)
ggarrange(g2, g3)
dev.off()
